Calibration of CEBPb and CEBPa - 02052020-AR

1 Introduction

In order to obtain the stoichiometry of CEBPa and CEBPb and how their relative abundances relate to PPARg abundance, cells were infected with a calibration construct. A linear model was fit to the calibration construct in order to model the 1:1 stoichiometry between CEBPa and CEBPb. Additionally, multiple over expression constructs were used to titrate the expression of CEBPa and CEBPb using dox.

Looking at your data, here are the 2 experiments I think would be needed. i think they should be relatively easy to do: 1) over-express CEBPB for 48 hours while knocking down CEBPA, plot calibrated PPARG vs CEBPB, if possible gate for low CEBPA expression; AND 2) Differentiate cells for 96 hours with DMI, induce CEBPB for 48 hours or control induction or CEBPB knockdown, measure calibrated CEBPA, CEBPB and PPARG, plot PPARG versus CEBPA for low, intermediate and high CEBPB, determine effect of equimolar CEBPB on PPARG expression.

With the first experiment, the goal is to test whether it is true that CEBPB cannot activate PPARG as well as CEBPA, but the complication is that we are always producing CEBPA when we OX-CEBPB. We have a lot of evidence for this being true that CEBPB cannot activate PPARG as well as CEBPA, but we need to cleanly show it.

With the second experiment, you have already shown that knocking down CEBPB late increases PPARG expression and CEBPA has the opposite effect. Here we just want to just show it again but this time quantify the levels of CEBPB and CEBPA.

 

In this experiment Dapi is channel 1, PPARg is channel 2, CEBPB is channel 3, and CEBPA is channel 4.

1.1 Import libraries

# Import libraries
library(tidyverse)
library(viridis)
library(ggridges)
library(ggpointdensity)
library(cowplot)

1.2 Initial Data Import

#Import MACKtrack import function from Fluffy
source('/Volumes/labdata/kyle/R functions/MACKimport.R')

exp_1 <- MACKimport('/Volumes/labdata2/atefeh/02052020-AR/Tracking/AllData_R.mat')
exp_2 <- MACKimport('/Volumes/labdata2/atefeh/02052020-AR/tracking2/AllData_R.mat')

saveRDS(exp_1, "exp_1.RDS")
saveRDS(exp_2, "exp_2.RDS")

Channels:

1. Dapi
2. Pparg (YFP)
3. Cebpb (mCherry)
4. Cebpa (Cy5)

2 Experiment 1

With the first experiment, the goal is to test whether it is true that CEBPB cannot activate PPARG to the same degree that CEBPA can. The complication is that we are always producing CEBPA when we OX-CEBPB. We have a lot of evidence for this being true that CEBPB cannot activate PPARG as well as CEBPA, but we need to cleanly show it. To overcome this complication, we over-expressed CEBPB for 48 hours while knocking down CEBPA.

In the analysis below we will show: * calibrated PPARG vs CEBPB, if possible gate for low CEBPA expression

2.1 Import and clean data

exp_1 <- readRDS("exp_1.RDS")
exp_1 <- exp_1 %>% mutate(Condition = gsub('0x2E','.', Condition),
                          Condition = gsub("C_", "C_NA_", Condition)) %>%
  separate(Condition, into = c("cell","siRNA","dox","rep"), sep = "_", remove = FALSE)

2.2 Filter out very low and very high integrated Dapi cells

ggplot(exp_1, aes(x = log2(MeanNuc1), y = Condition)) + geom_density_ridges2(scale = 5) +
  xlim(5,12) + theme_bw() +
  geom_vline(xintercept = c(6,10.5), color = "red", linetype = "dashed")

exp_1 <- exp_1 %>% filter(log2(MeanNuc1) > 6, log2(MeanNuc1) < 10.5)

2.3 Density plots for CEBPa, CEBPb, and PPARg for all 4 plates

2.4 Fitting linear model for 1:1 stoichiometry of CEBPa and CEBPb

# Calibration region
p1 <- exp_1 %>% filter(cell == "C") %>%
  #sample_n(1E5) %>%
  ggplot(., aes(x = log2(MeanNuc3), 
                y = log2(MeanNuc4))) + 
  annotate('rect', 
           xmin = 11, 
           xmax = Inf, 
           ymin = 11, 
           ymax = Inf, 
           alpha = 0.2, 
           fill = 'red') +
  geom_point(alpha = 0.1, size = 0.25) +
  theme_bw() +
  xlim(5,NA) + ylim(5,NA) +
  xlab('log2(Cebpb)') + ylab('log2(Cebpa)')

# Fit linear model
calibration <- exp_1 %>% filter(cell == "C",
                                log2(MeanNuc3) > 11,
                                log2(MeanNuc4) > 11) %>%
  lm(log2(MeanNuc3) ~ log2(MeanNuc4), data = .)

# Calibration region with linear model
p2 <- exp_1 %>% filter(cell == "C",
                 log2(MeanNuc3) > 11,
                 log2(MeanNuc4) > 11) %>%
  ggplot(., aes(x = log2(MeanNuc3), y = log2(MeanNuc4))) + geom_point() +
  geom_abline(intercept = calibration$coefficients[1],
              slope = calibration$coefficients[2],
              color = 'red', size = 3, alpha = 0.5) +
  theme_bw() + xlim(NA,15) + ylim(NA,15)

plot_grid(p1,p2, labels = 'AUTO')

2.5 Pairwise plots of CEBPa and CEBPb with linear model

exp_1 %>% filter(cell != 'C') %>% 
  rename(Pparg = MeanNuc2, 
         Cebpb = MeanNuc3, 
         Cebpa = MeanNuc4) %>%
  ggplot(., aes(x = log2(Cebpa), 
                y = log2(Cebpb), 
                colour = log2(Pparg))) +
  geom_point(size = 0.2, alpha = 0.2) + 
  facet_wrap(siRNA~dox) + 
  geom_abline(intercept = calibration$coefficients[1], 
              slope = calibration$coefficients[2],
              color = 'red',
              size = 1,
              linetype = 'dashed') +
  scale_color_viridis(direction = -1) +
  ylim(2,14) + xlim(2, 14) +
  theme_bw() +
  theme(text = element_text(size=10),
          legend.position="bottom")

Here CEBPa and CEBPb are plotted together with the linear model of the 1:1 stoichiometry between the two proteins from the calibration construct. PPARg abundance is represented by the color gradient, with darker blue representing high PPARg.

exp_1 %>% filter(cell != 'C') %>% 
  rename(Pparg = MeanNuc2, 
         Cebpb = MeanNuc3, 
         Cebpa = MeanNuc4) %>%
  ggplot(., aes(x = log2(Cebpa), 
                y = log2(Cebpb), 
                colour = log2(Pparg))) +
  geom_bin2d(bins = 100) + 
  facet_wrap(siRNA~dox) + 
  geom_abline(intercept = calibration$coefficients[1], 
              slope = calibration$coefficients[2],
              color = 'red',
              size = 1,
              linetype = 'dashed') +
  scale_fill_viridis(direction = 1) +
  ylim(2,14) + xlim(2, 14) +
  theme_bw() +
  theme(text = element_text(size=10),
          legend.position="bottom")

Density plot showing where most of the cell density is from the previous plot.

p1 <- exp_1 %>% 
  rename(Pparg = MeanNuc2, 
         Cebpb = MeanNuc3, 
         Cebpa = MeanNuc4) %>%
  filter(log2(Cebpa) > 4,
         log2(Cebpa) < 11,
         log2(Cebpb) > 5) %>%
  mutate(interval = cut_interval(log2(Cebpa), n = 5)) %>%
  ggplot(., aes(x = log2(Cebpa), 
                fill = interval)) +
  geom_histogram(bins = 500) + theme_bw() +
  theme(text = element_text(size=10),
          legend.position="top") +
  xlim(4,11)

p2 <- exp_1 %>% 
  rename(Pparg = MeanNuc2, 
         Cebpb = MeanNuc3, 
         Cebpa = MeanNuc4) %>%
  filter(log2(Cebpa) > 4,
         log2(Cebpa) < 11,
         log2(Cebpb) > 5) %>%
  mutate(interval = cut_interval(log2(Cebpa), n = 5)) %>%
  filter(cell != 'C',
         siRNA != 'NCsiRNA') %>% 
  ggplot(., aes(x = Cebpb,
                y = Pparg)) +
  geom_point(size = 0.2, 
             alpha = 0.2,
             aes(color = interval)) + 
  facet_grid(dox~interval) +
  ylim(NA,1000) + xlim(0, 15000) +
  theme_bw() +
  theme(text = element_text(size=10),
          legend.position="none")

p3 <- exp_1 %>% 
  rename(Pparg = MeanNuc2, 
         Cebpb = MeanNuc3, 
         Cebpa = MeanNuc4) %>%
  filter(log2(Cebpa) > 4,
         log2(Cebpa) < 11,
         log2(Cebpb) > 5) %>%
  mutate(interval = cut_interval(log2(Cebpa), n = 5)) %>%
  filter(cell != 'C',
         siRNA == 'NCsiRNA') %>% 
  ggplot(., aes(x = Cebpb,
                y = Pparg)) +
  geom_point(size = 0.2, 
             alpha = 0.2,
             aes(color = interval)) + 
  facet_grid(dox~interval) +
  ylim(NA,1000) + xlim(0, 15000) +
  theme_bw() +
  theme(text = element_text(size=10),
          legend.position="none")

plot_grid(p1,p2,p3, 
          ncol = 1, 
          labels = "AUTO", 
          rel_heights = c(3,8,2),
          axis = 'l',
          align = 'v'
          )

A. Histogram with 5 intervals of Cebpa abundance that is used in the following figures.

B. Pairwise plots between Cebpb and Pparg from cells that were treated with siRNA targeting Cebpa. Plots are separated based on Cebpa expression (A) and dox concentration.

C. Same scheme as B, but these cells were treated with non-targeting siRNA.

2.6 Mean Pparg abundance across intervals of Cebpa and Cebpb abundances

p1 <- exp_1 %>% 
  rename(Pparg = MeanNuc2, 
         Cebpb = MeanNuc3, 
         Cebpa = MeanNuc4) %>%
  filter(log2(Cebpa) > 4,
         log2(Cebpa) < 11,
         log2(Cebpb) > 5,
         log2(Pparg) > 4,
         log2(Pparg) < 10,
         cell != 'C',
         siRNA != 'NCsiRNA') %>%
  mutate(cebpa_interval = cut_interval(log2(Cebpa), n = 5),
         cebpb_interval = cut_interval(log2(Cebpb), n = 50)) %>%
  unite(col = 'condition', siRNA, dox, sep = '_') %>%
  group_by(cebpa_interval, cebpb_interval) %>%
  summarise(mean_pparg = mean(Pparg, na.rm = T),
            n = n(),
            sd = sd(Pparg, na.rm = T),
            sem = sd / sqrt(n)) %>%
  filter(n >= 10) %>%
  ggplot(., aes(x = cebpb_interval,
                y = mean_pparg,
                color = cebpa_interval,
                group = cebpa_interval)) +
  geom_path() + 
  geom_point() +
  geom_errorbar(aes(ymin = mean_pparg - sem, 
                    ymax = mean_pparg + sem),
                width = 0) +
  ylim(0,500) +
  theme_classic() +
  theme(text = element_text(size=10),
          legend.position="none",
        axis.text.x = element_text(angle = 45, hjust = 1, size = 5))



p2 <- exp_1 %>% 
  rename(Pparg = MeanNuc2, 
         Cebpb = MeanNuc3, 
         Cebpa = MeanNuc4) %>%
  filter(log2(Cebpa) > 4,
         log2(Cebpa) < 11,
         log2(Cebpb) > 5,
         log2(Pparg) > 4,
         log2(Pparg) < 10,
         cell != 'C',
         siRNA != 'NCsiRNA') %>%
  mutate(`Cebpb/Cebpa` = Cebpb/Cebpa,
         cebpa_interval = cut_interval(log2(Cebpa), n = 5),
         ratio_interval = cut_interval(`Cebpb/Cebpa`, n = 50)) %>%
  unite(col = 'condition', siRNA, dox, sep = '_') %>%
  group_by(cebpa_interval, ratio_interval) %>%
  summarise(mean_pparg = mean(Pparg, na.rm = T),
            n = n(),
            sd = sd(Pparg, na.rm = T),
            sem = sd / sqrt(n)) %>%
  filter(n >= 10) %>%
  ggplot(., aes(x = ratio_interval,
                y = mean_pparg,
                color = cebpa_interval,
                group = cebpa_interval)) +
  geom_point() +
  geom_vline(xintercept = 14, linetype = 'dashed') +
  geom_path() + 
  geom_errorbar(aes(ymin = mean_pparg - sem, 
                    ymax = mean_pparg + sem),
                width = 0) +
  ylim(0,500) +
  theme_classic() +
  theme(text = element_text(size=10),
          legend.position="none",
        axis.text.x = element_text(angle = 45, hjust = 1, size = 5))

legend <- get_legend(
  exp_1 %>% 
  rename(Pparg = MeanNuc2, 
         Cebpb = MeanNuc3, 
         Cebpa = MeanNuc4) %>%
  filter(log2(Cebpa) > 4,
         log2(Cebpa) < 11,
         log2(Cebpb) > 5,
         log2(Pparg) > 4,
         log2(Pparg) < 10,
         cell != 'C',
         siRNA != 'NCsiRNA') %>%
  mutate(cebpa_interval = cut_interval(log2(Cebpa), n = 5)) %>%
    ggplot(., aes(x = Pparg, fill = cebpa_interval)) +
    geom_histogram() +
    theme(legend.position = 'bottom')
)

plots <- plot_grid(p1,p2,
                   axis = 'b',
                   align = 'h',
                   labels = 'AUTO')

plot_grid(plots,legend,
          ncol = 1,
          rel_heights = c(2,0.4))

1. Pairwise plot of mean Pparg abundance vs segments of Cebpb abundance (50 segments of equal step sizes, similar to A in previous plot) color coded by Cebpa abundance intervals (see A in previous plot). Error bars are SEM.

2. Similar to A except that the x-axis is an interval is of log2(Cebpb / Cebpa). Interestingly, when the ratio of Cebpb / Cebpa is around 15, mean Pparg levels begin to drop (dashed line). In fact, for the highest levels of Cebpa (purple) the ratio never crosses past this point.

3 Experiment 2

3.1 Import and clean data

#cell_media_siRNA_dox_rep

exp_2 <- readRDS("exp_2.RDS")
exp_2 <- exp_2 %>% mutate(
  Condition = gsub('0x2E', '.', Condition),
  Condition = gsub("constant_d", "constant_NA_d", Condition),
  Condition = gsub("refreshed_d", "refreshed_NA_d", Condition),
  Condition = gsub("cebpbsiRNA_", "cebpbsiRNA_dox0_", Condition),
  Condition = gsub("C_R", "C_NA_NA_NA_R", Condition),
  Condition = gsub("C_d", "C_NA_NA_d", Condition)
) %>%
  separate(
    Condition,
    into = c("cell", "media", "siRNA", "dox", "rep"),
    sep = "_",
    remove = FALSE
  ) 

3.2 Filter out very low and very high integrated Dapi cells

ggplot(exp_2, aes(x = log2(MeanNuc1), y = Condition)) + geom_density_ridges2(scale = 5) +
  xlim(5,12) + theme_bw() +
  geom_vline(xintercept = c(7,10.5), color = "red", linetype = "dashed")

exp_2 <- exp_2 %>% filter(log2(MeanNuc1) > 7, log2(MeanNuc1) < 10.5)

3.3 Density plots for CEBPa, CEBPb, and PPARg

3.4 Fitting linear model for 1:1 stoichiometry of CEBPa and CEBPb

# Calibration region
p1 <- exp_2 %>% filter(cell == "C",
                       dox != "NA") %>%
  #sample_n(1E5) %>%
  ggplot(., aes(x = log2(MeanNuc3), 
                y = log2(MeanNuc4))) + 
  annotate('rect', 
           xmin = 11, 
           xmax = Inf, 
           ymin = 11, 
           ymax = Inf, 
           alpha = 0.2, 
           fill = 'red') +
  geom_point(alpha = 0.1, size = 0.25) +
  theme_bw() +
  xlim(5,NA) + ylim(5,NA) +
  xlab('log2(Cebpb)') + ylab('log2(Cebpa)')

# Fit linear model
calibration <- exp_1 %>% filter(cell == "C",
                                dox != "NA",
                                log2(MeanNuc3) > 11,
                                log2(MeanNuc4) > 11) %>%
  lm(log2(MeanNuc3) ~ log2(MeanNuc4), data = .)

# Calibration region with linear model
p2 <- exp_2 %>% filter(cell == "C",
                       dox != "NA",
                 log2(MeanNuc3) > 11,
                 log2(MeanNuc4) > 11) %>%
  ggplot(., aes(x = log2(MeanNuc3), y = log2(MeanNuc4))) + geom_point() +
  geom_abline(intercept = calibration$coefficients[1],
              slope = calibration$coefficients[2],
              color = 'red', size = 3, alpha = 0.5) +
  theme_bw() + xlim(NA,15) + ylim(NA,15)

plot_grid(p1,p2, labels = 'AUTO')

3.5 Pairwise plots of CEBPa and CEBPb with linear model

exp_2 %>% filter(cell != 'C') %>% 
  rename(Pparg = MeanNuc2, 
         Cebpb = MeanNuc3, 
         Cebpa = MeanNuc4) %>%
  ggplot(., aes(x = log2(Cebpa), 
                y = log2(Cebpb), 
                colour = log2(Pparg))) +
  geom_point(size = 0.2, alpha = 0.2) + 
  facet_wrap(siRNA~dox) + 
  geom_abline(intercept = calibration$coefficients[1], 
              slope = calibration$coefficients[2],
              color = 'red',
              size = 1,
              linetype = 'dashed') +
  scale_color_viridis(direction = -1) +
  ylim(2,16) + xlim(2, 16) +
  theme_bw() +
  theme(text = element_text(size=10),
          legend.position="bottom")

Here CEBPa and CEBPb are plotted together with the linear model of the 1:1 stoichiometry between the two proteins from the calibration construct. PPARg abundance is represented by the color gradient, with darker blue representing high PPARg.

exp_2 %>% filter(cell != 'C') %>% 
  rename(Pparg = MeanNuc2, 
         Cebpb = MeanNuc3, 
         Cebpa = MeanNuc4) %>%
  ggplot(., aes(x = log2(Cebpa), 
                y = log2(Cebpb), 
                colour = log2(Pparg))) +
  geom_bin2d(bins = 100) + 
  facet_wrap(siRNA~dox) + 
  geom_abline(intercept = calibration$coefficients[1], 
              slope = calibration$coefficients[2],
              color = 'red',
              size = 1,
              linetype = 'dashed') +
  scale_fill_viridis(direction = 1) +
  ylim(2,16) + xlim(2, 16) +
  theme_bw() +
  theme(text = element_text(size=10),
          legend.position="bottom")

Density plot showing where most of the cell density is from the previous plot.

p1 <- exp_2 %>% 
  rename(Pparg = MeanNuc2, 
         Cebpb = MeanNuc3, 
         Cebpa = MeanNuc4) %>%
  filter(log2(Cebpa) > 4,
         log2(Cebpa) < 12) %>%
  mutate(interval = cut_interval(log2(Cebpa), n = 5)) %>%
  ggplot(., aes(x = log2(Cebpa), 
                fill = interval)) +
  geom_histogram(bins = 500) + theme_bw() +
  theme(text = element_text(size=10),
          legend.position="top") +
  xlim(4,12)

p2 <- exp_2 %>% 
  rename(Pparg = MeanNuc2, 
         Cebpb = MeanNuc3, 
         Cebpa = MeanNuc4) %>%
  filter(log2(Cebpa) > 4,
         log2(Cebpa) < 12,
         log2(Cebpb) > 4) %>%
  mutate(interval = cut_interval(log2(Cebpa), n = 5)) %>%
  filter(cell != 'C',
         siRNA != 'NCsiRNA') %>% 
  ggplot(., aes(x = Cebpb,
                y = Pparg)) +
  geom_point(size = 0.2, 
             alpha = 0.2,
             aes(color = interval)) + 
  facet_grid(dox~interval) +
  ylim(NA,1000) + xlim(0, 15000) +
  theme_bw() +
  theme(text = element_text(size=10),
          legend.position="none")

p3 <- exp_2 %>% 
  rename(Pparg = MeanNuc2, 
         Cebpb = MeanNuc3, 
         Cebpa = MeanNuc4) %>%
  filter(log2(Cebpa) > 4,
         log2(Cebpa) < 12,
         log2(Cebpb) > 4) %>%
  mutate(interval = cut_interval(log2(Cebpa), n = 5)) %>%
  filter(cell != 'C',
         siRNA == 'NCsiRNA') %>% 
  ggplot(., aes(x = Cebpb,
                y = Pparg)) +
  geom_point(size = 0.2, 
             alpha = 0.2,
             aes(color = interval)) + 
  facet_grid(dox~interval) +
  ylim(NA,1000) + xlim(0, 15000) +
  theme_bw() +
  theme(text = element_text(size=10),
          legend.position="none")

plot_grid(p1,p2,p3, 
          ncol = 1, 
          labels = "AUTO", 
          rel_heights = c(3,8,2),
          axis = 'l',
          align = 'v'
          )

A. Histogram with 5 intervals of Cebpa abundance that is used in the following figures.

B. Pairwise plots between Cebpb and Pparg from cells that were treated with siRNA targeting Cebpa. Plots are separated based on Cebpa expression (A) and dox concentration.

C. Same scheme as B, but these cells were treated with non-targeting siRNA.

3.6 Mean Pparg abundance across intervals of Cebpa and Cebpb abundances

p1 <- exp_2 %>% 
  rename(Pparg = MeanNuc2, 
         Cebpb = MeanNuc3, 
         Cebpa = MeanNuc4) %>%
  filter(log2(Cebpa) > 4,
         log2(Cebpa) < 12,
         !is.na(log2(Cebpb)),
         cell != 'C',
         media == 'constant',
         siRNA == 'NA') %>%
  mutate(cebpa_interval = cut_interval(log2(Cebpa), n = 5),
         cebpb_interval = cut_interval(log2(Cebpb), n = 50)) %>%
  group_by(cebpa_interval, cebpb_interval) %>%
  summarise(mean_pparg = mean(Pparg, na.rm = T),
            n = n(),
            sd = sd(Pparg, na.rm = T),
            sem = sd / sqrt(n)) %>%
  filter(n >= 20) %>%
  ggplot(., aes(x = cebpb_interval,
                y = mean_pparg,
                color = cebpa_interval,
                group = cebpa_interval)) +
  geom_path() + 
  geom_point() +
  geom_errorbar(aes(ymin = mean_pparg - sem, 
                    ymax = mean_pparg + sem),
                width = 0) +
  ylim(0,NA) +
  theme_classic() +
  theme(text = element_text(size=10),
          legend.position="none",
        axis.text.x = element_text(angle = 45, hjust = 1, size = 5))



p2 <- exp_2 %>% 
  rename(Pparg = MeanNuc2, 
         Cebpb = MeanNuc3, 
         Cebpa = MeanNuc4) %>%
  filter(log2(Cebpa) > 4,
         log2(Cebpa) < 12,
         cell != 'C',
         media == 'constant',
         siRNA == 'NA') %>%
  mutate(`Cebpb/Cebpa` = Cebpb/Cebpa,
         cebpa_interval = cut_interval(log2(Cebpa), n = 5),
         ratio_interval = cut_interval(`Cebpb/Cebpa`, n = 50)) %>%
  group_by(cebpa_interval, 
           ratio_interval) %>%
  summarise(mean_pparg = mean(Pparg, na.rm = T),
            n = n(),
            sd = sd(Pparg, na.rm = T),
            sem = sd / sqrt(n)) %>%
  filter(n >= 20) %>%
  ggplot(., aes(x = ratio_interval,
                y = mean_pparg,
                color = cebpa_interval,
                group = cebpa_interval)) +
  geom_point() +
  geom_vline(xintercept = 15, linetype = 'dashed') +
  geom_path() + 
  geom_errorbar(aes(ymin = mean_pparg - sem, 
                    ymax = mean_pparg + sem),
                width = 0) +
  ylim(0,NA) +
  theme_classic() +
  theme(text = element_text(size=10),
          legend.position="none",
        axis.text.x = element_text(angle = 45, hjust = 1, size = 5))


legend <- get_legend(
  exp_2 %>% 
  rename(Pparg = MeanNuc2, 
         Cebpb = MeanNuc3, 
         Cebpa = MeanNuc4) %>%
  filter(log2(Cebpa) > 4,
         log2(Cebpa) < 12,
         log2(Cebpb) > 4,
         log2(Pparg) > 4,
         log2(Pparg) < 12,
         cell != 'C',
         siRNA != 'NCsiRNA') %>%
  mutate(cebpa_interval = cut_interval(log2(Cebpa), n = 5)) %>%
    ggplot(., aes(x = Pparg, fill = cebpa_interval)) +
    geom_histogram() +
    theme(legend.position = 'bottom')
)

plots <- plot_grid(p1,p2,
                   axis = 'b',
                   align = 'h',
                   labels = 'AUTO')

plot_grid(plots,legend,
          ncol = 1,
          rel_heights = c(2,0.4))

1. Pairwise plot of mean Pparg abundance vs segments of Cebpb abundance (50 segments of equal step sizes, similar to A in previous plot) color coded by Cebpa abundance intervals (see A in previous plot). Error bars are SEM.

2. Similar to A except that the x-axis is an interval is of log2(Cebpb / Cebpa). Interestingly, when the ratio of Cebpb / Cebpa is around 15, mean Pparg levels begin to flatten out, same ratio as in the previous experiment!!